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^ : ABSTRACT 

In this paper we derive equations describing dynamics and stationary con- 
figurations of a twisted fully relativistic thin accretion disc around a slowly 
0^ ' rotating black hole. We assume that the inclination angle of the disc is small 

m ; 

i and that the standard relativistic generalisation of the a model of accre- 

. tion discs is valid when the disc is flat. We find that similarly to the case 

' of non-relativistic twisted discs the disc dynamics and stationary shapes can 

be determined by a pair of equations formulated for two complex variables 
describing orientation of the disc rings and velocity perturbations induced by 
. the twist. 

. . .' We analyse analytically and numerically the shapes of stationary twisted 

configurations of accretion discs having non-zero inclinations with respect to 
the black hole equatorial plane at large distances r from the black hole. It 
is shown that the stationary configurations depend on two parameters - the 
viscosity parameter a and the parameter 6 =^ where is the opening 

angle ((5* ~ h/r, where h is the disc half thickness and r is large) of a flat disc 
and a is the black hole rotational parameter. When a > and (5^1 the shapes 
depend drastically on value of a. When a is small the disc inclination angle 
oscillates with radius with amplitude and radial frequency of the oscillations 
dramatically increasing towards the last stable orbit, Rms- When a has a 
moderately small value the oscillations do not take place but the disc does 
not align with the equatorial plane at small radii. The disc inclination angle is 
either increasing towards Rms or exhibits a non-monotonic dependence on the 
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radial coordinate. Finally, when a is sufficiently large the disc aligns with the 
equatorial plane at small radii. When a < the disc aligns with the equatorial 
plane for all values of a. 

The results reported here may have implications for determining structure 
and variability of accretion discs close to Rms as well as for modelling of 
emission spectra coming from different sources, which are supposed to contain 
black holes. 

Key words: accretion, accretion discs; hydrodynamics; black hole physics; 
relativity; binaries: close; galaxies: nuclei; celestial mechanics 



1 INTRODUCTION 



Twisted accretion discs around rotating black holes are used for explanations of differ- 
ent observational phenomena. For example they a re invoked to explain jet precession ob- 



served in certain active galactic nuclei (AGN) (e.g. ICaproni. Mosquera Cuesta Sz Abraham 



(1200. 



(e.g. 



) and references therein ) , warping of the accretion disc in the maser gala xy NGC 4258 



Papaloizou. Terquem fc LinI (Il998[ ): 



Caproni et al 



tures of emission line pro files c oming from AGN 



Wu. Chen &: YuanI ( 120101 )) etc.. iFerreira &: Ogilvid fcOOSl . |2009[ ) proposed that shear veloci- 



fcOOTh: 



Bachev 



1999f l: 



Martini (|2008[)'). certain fea - 

(bood); 



Cadez et al. 



ties induced by the disc twist (see below) provide a mechanism of excitation of high frequency 
quasi-periodic oscillations observed in many astronomical objects. Light curves, forms of 
emission lines of precessing twisted tori around a Kerr black hole as well as their shapes as 
seen from large distances have recently been investigated numerically by iDexter fc Fragile 



(1201 ih 



Theoretical studies of thin twisted discs started from the seminal paper 



Bardeen fc Petterson 



( 1l975l ). who proposed a simple model of a stationary twisted disc, where a twisted disc was 
considered as a collection of circular rings interacting with each other due to viscous forces. 
In frameworks of this model Bardeen and Petterson came to the conclusion that a station- 
ary twisted disc inclined with respect to the black hole equatorial plane at large distances 
from the black hole aligns with this plane at smaller radii - the effect called later as "the 
Bardeen- Petterson effect". This alignment was explained as being due to the presence of 
'gravitomagnetic force' acting from the side of a rotating black hole on the disc rings. This 
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force causes precession of orbital planes of free particles around a direction of the black hole 
spin and breaks the spherical symmetry of the problem. The Bardeen and Pe tterson approach 



was 



19781); 



a ter developed and generalised on non-s tationary twisted discs by e.g. IPettersonI (119771 . 



Hatchett. Begelman fc SarazinI ( 1l98ll ). 



In frameworks of hydrodynamical theory of perturbations iPapaloizou &: Pringld ( 119831 ) 
pointed out that the simple model developed in the previous works is, in fact, inconsistent 
since it does not conserve all components of the angular momentum content of an accretion 
disc. They showed that a self-consistent model must, necessarily, contain perturbations of 
the disc density and velocity fields induced by the disc twist, which determine de v iation s of 
the gas particles motion in the disc from the circular one (see llvanov &: lUarionovl (Il997l ) for 
a qualitative explanation of this effect) . In particular, the velocity perturbations should be 
odd with respect to t h e disc vertical coordinate thus having the structure of a shear flow. 



Papaloizou &: Pringld ( 119831 ) also found that when a sufficiently viscous disc is considered 



in the classical Newtonian gravitational field of a point mass a typical alignment scale of a 



stationary twisted configur ation decreases wi t h dec rease of the 



Shakura &: SunyaevI (1l973l ) 



a parameter (see also e.g. iKumar fc Pringld (119851 )) and that non-stationary propagation 
of twisted disturbances through the disc has a diffus i ve cha racter with a characteristic time 
scale being also proportional to a (also e.g. iKumarl ( Il990l )) contrary to what was claimed 
in the previous studies. These two effects are determined by the well known degeneracy of 
the Keplerian poten tial, which has orb i tal an d epicyclic frequencies equal to each other. It 



was later shown by 



Papaloizou fc LinI ( 1l995l ) that in the opposite case of a low viscosity 



Keplerian disc propagation of non-stationary twisted disturbances through the disc has a 
wave-like character with a typical speed of the order of speed of sound in the disc. These 
studies considered the gravitational field of a black hole as a field of a Newtonian point mass 
with the gravitomagnetic force being treated as a classical force causing precession of the 
dis c rings. 



Ivanov &: Illarionovl ( 119971 ) considered post-Newtonian corrections to equations describ- 
ing velocity perturbations. They found that when a low viscosity stationary twisted disc 
is considered and the disc gas is rotating in the same sense as a black hole the Bardeen- 
Petterson effect does not take place. Instead, there are radial oscillations of the disc inclina- 
tion angle with amplitude increasing with decrease of the distanc e from the black hole. This 



effect w as later confirmed by time dependent calculations made by 



Lubow. Ogilvie &: Pringle 



( I2OO2I ). llvanov fc Illarionovl (119971 ) also provided a qualitative explanation of this effect and 
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a general criterion of appearance of these oscillations. Namely, they take place in a low vis- 
cosity stationary twisted disc when the signs of the apsidal and nodal precessions are the 
same. In the opposite case the disc aligns with the equatorial plane of the black hole for any 
value of viscosity. For example, the latter situation is realised when the black hole rotates 
in the direction opposite to the direction of the orbital motion of the gas in the disc, see 



Ivanov Sz lUarionovl ( 119971 ) and below or when a Ne wtonian twisted disc around a r nassiv e 



Ivanov. Papaloizou &: PolnarevI (Il999[ ). 



binary star or a binary black hole is considered, e.g. 
We are going to show below that this criterion is applicable to our fully relativistic problem 
as well and that whether the disc aligns with the equatorial plane or not depends on the 
direction of rotation of the black hole relative to the orbital motion. 
Additionally, 



Ivanov fc lUarionovl f|l997n: 



linear analysis made by 



Demianski fc Ivanovl f|l997l) 



Papaloizou fc Pringld (119831 ) and iPapaloizou fc Lid ( 119951 ) can be 



showed that the 



extended to the case of disc inclination angles, which are larger than the disc opening angle 
6^: ~ h/r ^ 1, where h is the disc half thickness and r is the radial coordinate, and the 
precise definition of the angle 5=,, is made in such way that it is a constant independent on 
r, see equations (l37l) and ( l38l) below. This was done by employi n g a formalis m based on the 

(119771 . 119781 ). An analogous 



Petterson 



so-called "twisted coordinate" system intro duced b y 
coordinate system has later been used by Ogilvie ( 19991 ) to consider the disc inclination 
angles of order of unity and, accordingly, to construct a non-linear theory of evolution of 
sufficiently viscous twisted discs. 

As we described above, all previous formalisms of description of the thin twisted discs 



were based either on the classical or post-Newtonian treatments of the probleno. It is, 
however, important to consider the problem in full General Relativity. This is especially 
crucial either for the discs having a small value of a ^ 5* or in the case of slowly rotating 
black holes having their rotational parameter a ^ 1, where a stationary accretion disc can 
be twisted at scales comparable to the black hole gravitational radius, Vg. A formalism based 
on General Relativity can also be very useful for a quantitative modelling of emission spectra 
coming from twisted disc as well as for accurate studying of many other physical processes 
occurring in such discs like e.g. their self-irradiation. 



^ Note that there is a number of misprints in 



Demianski &: Ivanov! 11997). The final equations of this paper are, however, 



correct 
2 



Note, however, that dyn a mical models of fu l ly rel ativistic accretion tori have been recently considered by numerical means, 



e.g. 



Fragile fc Blaej l l200Sri : 



Dexter &: Fragile 



1 2OI1I') . 
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In this paper we construct a fully relativistic theory of thin twisted discs at radii larger 
than the radius of the last stable orbit around slowly rotating black holes, treating the effect 
of the black hole rotation in the linear approximation only. This essentially means that only 
one term in general equations of motion (proportional to the rotational parameter a and 
describing the gravitomagnetic force) can be taken into account and all quantities entering 
in our final equations describing dynamics of twisted discs can be defined with help of the 
Schwarzschild metric of a non-rotating black hole. We make the usual assumptions of small- 
ness of the disc opening angle 5*, the disc inclination angle (3 with respect to the equatorial 
plane and its logarithmic derivative, r-^(3 ^ 1. Additi onally, we neglec t self-g ravity of the 



PettersonI fll977l . 



19781 ) to the rela- 



discs. We generalise the twisted coordinate system of 
tivistic case and write down our general equations of motion in this generalised coordinate 
system. We show that these equations can be naturally split on two different subsets: 1) a set 
of equations describing the standard flat disc accretion and 2) a set of equations determining 
dynamics of variables relate d to the disc twist. As a s et of the standard equations we use 
equations of the well known iNovikov fc Thornd (119731) ac cretion disc model, which is the 



relativistic generalisation of the IShakura fc SunyaevI ( 119731 ) a-disc model. Equations of the 
set 2) are used to derive a pair of equations (l60l) and ( 161]) fully describing the dynamics and 
stationary configurations of twist ed discs under our approxirn ations. Similar to the case of 
Newtonian twisted discs (see e.g. iDemianski fc Ivanovl (119971 )) these equations can be for- 
mulated for two dependent complex variables - W(r, t) = /3e'*', where 7 is the second Euler 
angle determining longitude of ascending nodes of the disc rings, and B(r, t) determining 
the form of velocity perturbations in the disc. 

In the second part of the paper we study in detail shapes of stationary configurations, 
which are described by equation ( l62l) . We show that this equation contains only two in- 
dependent parameters - a assumed to be a constant throughout the disc and 6 = 6^:/y/a. 
When the disc gas rotates in the sarn e dire ction as the black hole (a > 0), for the gas 



pressure dominated 



Novikov fc Thornd (119731 ) models we find that there are three possible 



qualitatively different shapes of the twisted discs depending on values of a. When a is very 
small the disc exhibits oscillations of the inclination angle, which have a different charac- 
ter of amplitude and radial frequency behaviour at large distances from the black hole and 
close to the last stable orb it. The disc behaviour at large radii is the same as described by 



Ivanov fc lUarionovl (119971 ). On the other hand oscillations of the inclination angle close to 



the last stable orbit proceed in a different regime, with the amplitude and the radial fre- 
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quency strongly growing towards the last stable orbit. When a has moderately small values 
the oscillations are suppressed but the disc does not align with the black hole equatorial 
plane. The inclination angle is either monotonically growing towards the last stable orbit or 
exhibits a non-monotonic behaviour decreasing with decrease of r at large radii and increas- 
ing again towards the last stable orbit at smaller radii. Finally, when a is sufficiently large 
the Bardeen-Petterson effect takes place. 

In Section 2 we present our basic definitions, introduce the twisted coordinate system, 
describe fundamental equations of motions and the flat disc model used in our study. In 
Section 3 we derive our final equations ( l60l) and ( 16T1) describing the dynamics and stationary 
configuration of the twisted discs. Section 4 is devoted to analytical and numerical studies 
of the stationary twisted configurations. Note that Section 3 and Section 4 can be read 
independently. 

We use the natural systems of units throughout the paper setting the gravity constant 
and the speed of light to unity. The metric signature is (+, — , — , — ) and the usual summation 
rule over repeating indices is implied. 



2 BASIC DEFINITIONS AND EQUATIONS 
2.1 Metric 

As was mentioned in Introduction we assume that the black hole rotates slowly with the 
rotational parameter a <^ 1 and derive our dynamical equations taking into account only 
leading terms in a. Therefore, for our purposes, it suffices to consider the metric of a Kerr 
black hole in the Boyer-Lindquist coordinates taking into account only linear in a terms 

ds^ = {1- 2M/R)df - (1 - 2M/R)-^dR^ - R\d9^ + sin^Odcp^) + Aa—sin^O dcf) dt. (1) 

R 

It is easy to see that apart from the g^^ term the metric coincides wit h the metric of a non- 



rotati ng black hole written in the Schwarzschild coordinates, see e.g. 
Jl971f ). 



Zeldovich fc Novikov 



In the non-relativistic problems the dynamical equations for the twisted disc s take sim- 



PettersonI (11977 



19781 ). As 



plest form in the so-called "twisted coordinates" introduced by 
we discuss in Section 2.2 these coordinates are obtained from the cylindrical ones by rota- 
tion. In order to find a simple relativistic generalisation of the twisted coordinates we would 
like to bring the metric ([1]) into another, the so-called 'spatially isotropic' form, where the 
spatial part of the line element is proportional to the Cartesian line element. This is done 
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by the change of the radial coordinate 



R = Rj[l + 



M 
2R, 



(2) 



which brings ([T]) in the form 



1 - ^ 

2Ri 
1+2^ 



di' 



1 + 



2Ri 



{dU\ + R]de'' + Rjsin^^^^^ 



a sin^O 
+ 4 —dt d(p, 



Ri 1 + 



2Ri J 



(3) 



where the radial and time coordinates are expressed in units of M from now on. Introducing 
the cylindrical coordinates {r = Rj sin O^cp^z = Rj cos 0} we can rewrite ^ in an equivalent 
form 

ds^ = Kldf + 2ar^KiK3d(f)dt - K^{dr^ + dz^ + r^d(j)'^), (4) 

where 



1 - 



1 + — ' 

^ ^ 2Ri 



1 



1 



1 

2B.I 



2 ■ 



2Rj J i?j ^ 

The metric ([3]) induces an associated orthonormal tetrad of one-forms 

= Kidt + ar^Ksdcp, uj'' = K2dr, uj'*' = rK2d(f), uj^ = K2dz. 
We also use Cartesian coordinates {x, y, z) related to (r, 0, z) in the usual way. 



(5) 



(6) 



2.2 The twisted coordinate system 

The twisted coordinates (r, rf^,V',0 obtained from the Cartesian ones by rotation 
f T \ /lO o\/t\ 



\ e / 



C0S7 sm7 
—sin'y cos'^ (3 
\ Psiri'-^ —f3cos'y 1 J 



X 

y 

\ ^ I 



(7) 



where the Euler angles P{T,r) and 7('r, r) are functions to be determined. The inclination 
angle /3 is assumed to be small. It is evident that r^^ + = + z^. 

Since only rtw is used below we omit the index (tw) later on. It is also convenient to use 
9? = ■?/' + 7 in our expressions and introduce the variables 

"^1 = (3 cos 7, ^'2 = /3 sin 7 (8) 

instead of (3 and 7 and the quantities 

Z = /3sin7/' = ^isin(^-^2COS(^, U = Z, W = Z' , (9) 

where partial derivatives over r and r are denoted by dot and prime, respectively. 



8 V. V. Zhuravlev and P. B. Ivanov 

The orthonormal one-forms associated with the twisted coordinate system are given by 
the expressions 

= {Ki - ariK^d^U)dT + a^K^d^^Z - rW)dr + arK?Xr - iZ)d^ - arK^d^Zd^, (10) 



uj'' = -^K2UdT + ^2(1 - ^W)dr, 



(11) 



(12) 

u:^ = rK2UdT + rK2Wdr + Kid^, (13) 

see Appendix A for their relation to the coordinate forms ([6]). 

When /3 = /?' = 7' = the twisted coordinates are reduced to cyhndrical coordinates 
rotated with respect to the coordinates introduced above by angle 7 and the forms fllOtllSp 
coincide with the expressions ([6]) \ 

The adjoint basis vectors Sj are obtained from the duality condition u;*(ej) = 5* . Explic- 
itly, we have 

^ ^ rUdA , 



= (dr + ^Udr + -d^Ud^ 
Ki\ r 



= (-a^^d^Zdr + (1 + ^W)dr + -d^Wd^ - rWd, 
K2 \ Ki r 

-a^{r - iZ)dr - a^^rUdr + f - - a^^d^U] + ar^rUd^ 
Ki Ki \r Ki J ^ Ki ^ 



(14) 
(15) 
(16) 
(17) 



In the limit of large r the ex pressions ( ITOlfTTll tend to the corresponding expressions 



calculated for a flat space-time in 



Pettersonl fllQTSh 



We project all our dynamical quantities and equations of motion onto the basis (HUl - 
[T71) . In order to calculate covariant derivatives in this basis the connection coefficients, Tijk, 
should be used. Their explicit forms are given in Appendix A. 

For our purposes we need covariant derivative of a vector and covariant divergence of a 
tensor projected onto the bases fll0m7p : 

A^=e,(AO + n^A^ A^ = e,{A^^) + Tl^A'^+Ti^jV', (18) 

where (; ) stands for covariant derivative. Note that the indices of components of the tensor 



^ Note that when a = it is sufBcient to have /3' = 7' = for ||10I13P to coincide with JgJ. Clearly, this is due to the spherical 
symmetry of the Schwarzschild space-time. 
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quantities as well as the ones of the connection coefficients projected onto our orthonormal 
bases can be raised and lowered with help of the Minkowski metric rjik- 



2.3 Equations of motion 

Our equations of motion follow from the law of mass conservation 

ipU% = 0, (19) 

where p is the rest mass density, W are components of four velocity, and equality to zero of 
covariant derivative of the stress energy tensor, 

T-k = 0. (20) 

We consider the stress energy tensor of a viscous radiative fluid 

r'' = (e + p)U'U'' - pg''' + - Wq'' - U^q\ (21) 

where e and p are the energy density and pressure, respectively, g^^ and Tl^ are the compo- 
nents of the metric and viscosity tensors, are the components of radiation flux. We have 
e = p + eth, where the thermal energy eth is assumed to be much smaller than the rest energy 
p. The viscous stress tensor T^^ = 2770"*'^, where t] is dynamical viscosity. 



■ 1/ 

ik 



a 



-{U'.P^^ + U^.P^') - -W.P 

2 J 'J 3 d 



ik 



is the shear tensor and 



is the projection tensor. Note that 



ik 



-ik 



0, 



iTTk 



0. 



(22) 



(23) 



(24) 



Similar to the analysis of Newtonian twisted discs in the fully relativistic case in the 
linear approximation in angle f3 we can divide our equations of motion (IT^ and (I2U]) in 
two parts having different symmetries with respect to the coordinates C, and 1) a part 
describing an unperturbed flat disc and 2) a part describing the disc's perturbations asso- 
ciated with the disc's twist and warp. The ffist (background) part follows from standard 
models of relativistic flat disc accretion, where all dynamical quantities should be pro- 
jected onto the basis (11411171) with Z, U and W formally set to zero. In this paper we 
would like to consider the simplest possible case of a relativistic geometrically thin opti- 
cally thick stationary a— d i sc as a background mode l, and use, acc ordingly, the results o f 



Shakura fc SunyaevI (119731 ) . iNovikov fc Thornd (jl973[ ). hereafter NT. lPage fc Thornd (119741 ) 
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and 



Riffert fc Heroldl (119951 ). herea: 



in a way similar to what is done in 



ter RH, to describe it. The perturbation part is treated 



Demianski fc Ivanovl ( 119971 ). hereafter DI, it is discussed 



in Section 3 with some technical details relegated to Appendix B. 
2.4 The background model of a flat relativistic disc 

For our purposes it is sufficient to consider only a few basic properties of the fiat disc 
models. Additionally, since we consider the case of a slowly rotating black hole in this 
Section we neglect corrections due to a nonzero value of a and discuss only accretion discs 
around a Schwarzschild black hole, formally setting a = 0. However, certain basic expressions 
accounting for linear in a corrections are shown in Appendix B, for completeness. 

Equations describing the relativistic flat disc models can themselves be split in two parts: 
the ones responsible for disc's structure and behaviour determined by effects occurring on 
a dynamical time scale and the ones determining the disc's structure and behaviour related 
to a slow viscous time scale. 

Let us discuss the first group of equations temporarily neglecting viscous interaction 
between neighbouring rings of the disc. In this approximation the disc's gas is orbiting in 
the azimuthal direction with nearly geodesic four-velocity having only two components: If^ 
and U"^ = a/1 + {U'^y. We use, accordingly, the standard expressions for geodesic circular 
motion in the field of a Schwarzschild black hole 

U- = {R- 3)-^/^ =(^r-2 + ^y^\ U-= = K,Kl'\^''U^, (25) 

where we set R = K2{r)r using the fact that we consider the disc situated close to the 
equatorial plane of a black hole and assume that the disc's gas is orbiting in positive direction 
with respect to the angle if. An important quantity associated with the circular motion is 
its frequency with respect to a distant observer, Vt = From equation (|25|) we get 

n = = R~'/'. (26) 

The disc structure in the vertical direction is governed by equation of hydrostatic balance 



TV 



P \ 

where we neglect eth ^ p- 



2 



(27) 



Viscous interactions in the disc result in a slow drift of the disc's gas in the radial 
direction, and, accordingly, an additional radial component of four velocity. If", appears. 
Equations determining disc's properties related to this effect follow from the laws of conser- 
vation of the rest mass, energy and angular momentum, radiation transfer in the vertical 
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direction and properties of viscous interactions in the disc. They allow one to complete the 
set of equations of the disc structure and obtain explicit expressions for basic quantities 
describing the disc. Let us discuss some of these equations relevant for our purposes. 
The law of the rest mass conservation may be written in the form 



M 
2^ 



where the surface density 



T.KiKlrU'' 



dip, 



(28) 



(29) 



and the rate of flow of the rest mass through the disc, M, is a constant. 

In the fiat disc's models only (r y) and {rr)— conipone nts of the viscous stress tensor are 
important. As is shown in NT and lPage fc Thornd ( 1l974l ) (see also RH) the laws of energy 
and angular momentum conservation result in an explicit expression for integrated over the 
vertical coordinate (r(/9)-component of the stress tensor Q 



271 



where the bar stands hereafter for vertically integrated quantities, 

V6 73, {y-V3){3 + 2V2) 



D = 1 



— In 

y 2?/ 



(30) 



(31) 



and y = y R = \fK^. Note that D = at the marginally stable orbit, when R = R^s = 6, 
where the disc truncates. The value of is obtained from (El]): T*^^ = ^T"^^. From 
and ( l25l) it follows another expression for T'''^, 

f''' = lv§-m'U''/r. (32) 

Equating (02]) to ([30]) we get 

M ^-1/2 n 

(33) 



V= 



Taking into account that when r — )■ oo the expression in the brackets tends to unity we get 
the well known result that f? = |^ is a constant in the nonrelativistic limit. 

For our purposes we need an expression for kinematic viscosity v = rj/ p. \n this paper 
we employ the simplest possible a pproach assumi ng th at v does not depend on th e vertical 



coordinate ^ and use the standard 



Shakural ( 119721 ) and IShakura fc SunyaevI (119731 ) prescrip- 



■* Note that this component is defined in NT and RH with respect to another, nearly comoving, frame of basis vectors. The 
component in our frame can be obtained multiplying the NT and RH result by . Also, integration of this component in these 
papers is performed over the proper vertical length Improper = K2S,, see l(4]l. Therefore, in order to obtain a vertically integrated 
value of T^f the NT and RH result should be additionally divided by K2- 
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tion for an estimate of this quantity assuming that it is proportional to the product of a 
characteristic disc scale height and sound speed, Cs, in the disc 

U ~ aCshproper, (34) 

where hproper is the proper scale height of the disc related to a coordinate, h, of the disc upper 
boundary in our coordinate system as hproper = -^2^, and a is the usual Shakura-Sunyaev 
parameter. It is assumed to be a constant. From the hydrostatic balance equation f l27p it 
follows that Cg ~ \fPTp ~ U'^h/r. Substituting this relation to we finally define a in 
such a way that 

V = aK2U'^h^/r. (35) 

Using equations (l33ll and ( l35ll we obtain an important relation 

2 _ M / D \ 

~^\U^{U-fKlKr)- ^^^^ 
Finally, we are going to use the ratio 5{r) = h/r in our calculations. It can be easily 
obtained from the results of NT provided that one takes into account a correction to the 
vertical balance equation of NT found in RH. When the disc is gas pressure dominated and 
the main source of opacity is determined by the Thomson scattering, we have 

5(r) = 6,K',^'Kl^'\un-^/'D'/'r'/'', (37) 

and for the gas pressure dominated disc with the free-free processes giving the main source 
of opacity we get 

5{r) = 5,Kf^''Kl'\U^)-^"^''D^/^^r^/\ (38) 

where 5=,, ^ 1 is a constant. Note that in both cases 5 = at the marginally stable orbit. 
It is also important to note that the disc opening angle 6{r) so defined does not depend on 
whether we use the Boyer-Lindquist or the isotropic coordinates. Indeed, taking into account 
that hproper = -^^2^ and R = we have 6{r) = h/r = hproper I R- 

The expressions for the flat NT model listed above are strictly valid only when the radial 
drift velocity is small: |f/''"| ^ Cg. This inequality is, however, broken when x = R — Rms is 
small, and, when R is very close to the position of the last stable orbit |f/^| > Cg. Assuming 
that the NT expressions can also give order of magnitude estimates even when |f/''| ~ let 
us estimate a value of x, Xg, where we have \U^\ = Cg. 

The expression for \U^\ follows from equations ( l28l) and ( l36i) . 

|f/'-| = ^^iKfU^Un'VR, (39) 
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while the speed of sound may be estimated as Cg ~ If^S, see above. Assuming that the new 
radial variable x is small we have from ( 1311) 

D^^, (40) 

we represent the quantity 6 as 6 = Smsx"^^, where 6ms ~ f^*, and e = 1/5, 3/20 for the cases of 
the disc opacity dominated by the Thomson scattering and free-free processes, respectively, 
see (1371) and (l38l) . In the same limit we can set the values of all variables entering (139|) and 
the expression for Cg excepting 6 and D equal to their values at Rms- We get 

\U'\/cs ~ A,a6*x^^'-'\ A, = ^KIU^W^^/R ~ 2.3 ■ 10"^, (41) 

48 

where we use ( I40p . From this equation we obtain 



'his e quation tells that for typical values of a and S 

ms ~ 10 — 10 (see e.g. I Ivanov. Igumenshchev fc Novi 



(jl998l ). their equation (1) ) the value of Xg is quite small. 



3 DERIVATION OF EQUATIONS GOVERNING DYNAMICS AND 
STATIONARY SHAPE OF A TWISTED DISC 

3.1 Basic facts and simplifications 

As described above we use equations ( 12U]) written in the basis (I14m7p to obtain equations 
determining dynamics and stationary configurations of a twisted disc and we take into 
account that our assumption about smallness of the inclination angle (5 allows us to decouple 
the set of equation (!20|) onto two subsets: 1) the subset of equations describing the standard 
flat relativistic disc reviewed above (the "standard subset") and 2) an additional set of 
equations determining propagation of warped and twisted disturbances through the disc 
and shapes of non-planar stationary configurations (we call it later as a subset describing 
the disc twist). All equations of this additional subset depend harmonically on the angle ip 
(or if) and can be obtained from the general set (!20!) by a procedure employing symmetry of 
different terms in (|20|) with respect to the vertical coordinate ^ since in the approximation 
we use all terms in this subset are either even or odd with respect to the change ^ — ?■ —S,- 
Namely, let us call the equations of ( l20l) having the index i = r, r, y9 as "the horizontal part 
of equations" and the equation with i = as "the vertical equation". Then, all terms of the 
horizontal part having the odd symmetry and all terms of the vertical equation having the 
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even symmetry belong to the subset of equations describing the disc twist while all terms 
with the opposite symmetries belong to the standard subset. 

There are additional simplifications determined by the fact that we consider only the 
thin discs /i/r ^ 1 and slowly rotating black holes with a <^1. 

The first assumption allows us to consider only terms of order of h/r in the horizontal 
part and terms of zero and second order in h/r in the vertical equation. In fact, as we 
are going to show there is only one zero order term in the vertical equation describing 
precession of the disc rings due to the Lense-Thirring effect (the so-called Lense-Thirring or 
gravitomagnetic precession). It is, therefore, proportional to a. Thus, this zero order term 
may also be classified as a "small" one due to our second assumption. A further simplification 
follows from the fact that the pressure p, the thermal density eth and the dynamic viscosity 
rj are small in comparison with p: p,eth,ri ~ {h/rYp. Thus, these quantities enter in the 
horizontal part only in expressions containing differentiation over ^ in order to provide 
terms linear in h/r. In the vertical equation we set e = p due to the same assumption. This 
results in absence of eth in all our expressions used below, see Appendix B. 

The second assumption allows us to suppose that a characteristic evolution time scale of 
propagation of twisted and warped disturbances through the disc, tf^, is much larger than 
the disc dynamical time scale td ~ VL~^. Indeed, for a non-rotating black hole with a = the 
disc rings inclined with respect to each other can change their mutual orientations only due 
to interactions determined by either pressure or viscous forcej^. Due to our first assumption 
these interactions are small and the corresponding time scales must be proportional to an 
inverse power of h/r. For a slowly rotating black hole there is an additional characteristic time 
scale, tiT, determined by the gravitomagnetic precession of the rings. It may be estimated 
as tiT ~ a~^r^ ^ td, see e.g. equation (1B11I) . Due to the small ratios td/Uw and td/tir the 
time derivatives of the Euler angles, and, accordingly, the variable U defined in ([9]) can be 
neglected in all expressions entering in the horizontal p art. 



As was first noted by iPapaloizou fc Pringld ( 1l983l ). hereafter PP, for the Newtonian 



problem in order to get a self-consistent set of equations describing the disc twist one must 
take into account perturbations of velocity, gas density and pressure determined by the 
twisted disturbances in the disc. In our relativistic generalisation we must use perturbations 



^ Let us remind that we neglect self-gravity of the disc. 
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of four-velocity, and, accordingly, assume that 

W = Wq + v\ p = po + p,, p = Po+Pi (43) 

where the quantities with index (0) are determined by the equations of the flat disc model 
described above. Note that in a self-consistent approach we should consider perturbations of 
the dynamical viscosity, rj = rjQ + rji, and radiation flux, = qo + ql, as well. However, it can 
be shown that contribution of terms proportional to the unperturbed part of radiation flux, 
its perturbation and perturbation of the dynamical viscosity, is negligible in our equations at 
the required order of accuracy. This is related to the fact that the quantities rji and ^ are 
proportional to {h/r)"^ and cx {h/r)q^. As we discuss above the terms proportional to 

{h/rY give a contribution to the horizontal part of our equations being differentiated over ^ 
and in the vertical part being multiplied by terms of zero order in h/r. A direct examination 
of our set of the perturbed equations shows that this, in fact, does not happen in the linear 
approximation in angle /3. Also note that the pressure perturbation pi does not enter in our 
set of final equations describing a twisted disc, see the next Section. 

Since both f/* and Uq obey the same normalisation condition t/jf/* = UqiUq = 1 the 
velocity perturbations are orthogonal to Uq. v^Uq = 0, giving = {Uq/Uq)v'^. The 
"horizontal" velocity perturbations are odd functions of C, while the "vertical" 

perturbation is an even function of the same coordinate. 

The horizontal part of the velocity perturbation is directly determined from our equations 
for the disc twist, see below. Determination of the vertical part, req uires some additional 



consid eration. Namely, as was first shown for the Newtonian problem by lHatchett. Begelman fc Sarazin 
( 1981 ). consists of two parts: = d^/dt + f^, the first one is determined by a change 
of the coordinate ^ with time of a given gas element while the second one is due to the 
non-coordinate character of our basis vectors (11411171) . This results in evolution of projection 
of the velocity vector f* onto the vector even when a gas element has a fixed value of 
^ provided that either the Euler angles are changing with time or there is a non-zero drift 
velocity component Uq. In order to find vl let us assume that d^/dt = 0. In this case we 
have = oj^/ds, where lo^ is given by equation (1131) . and ds is the line element. Using this 
equation, taking into account that f/J" = Kidr/ds and Uq = K2dr/ds, see equations ( fTOjl 
and ( ITTl) . respectively, we get 

= rU^^U + rU^W. (44) 
In what follows we set d^/dt = 0. This requirement is related to a specific freedom of choice 
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between different twisted coordinate systems corresponding to the same physical situatior 
Indeed, in order to construct our twisted coordinate system we add two Euler angles to 
the set of the state variables describing our physical system while the number of dynamical 
equations remains the same. Therefore, in the twisted coordinate system there is a freedom 
of choice between different sets of dynamical variables satisfyin g the same equations, which 
must be fixed by an additional requirement. As was discussed in llvanov fc Papaloizoul (120081 ) 
the requirement d^/dt = Q fixes this freedom and leads to a choice of the most appropriate 
twisted coordinate system, where perturbations of the state variables are minimal!. 

Since the unperturbed parts of our variables and the perturbations enter in a different way 
in all equations below, we omit for simplicity the index (0) characterising the unperturbed 
variables from now on. 

The horizontal part of our equations and the vertical equation are shown in Appendix 
B, see equations f lBlUBSP and equation (]B4p . respectively. They have a rather complicated 
structure. However, they can be further simplified. Indeed, a simple analysis shows that all 
terms in these equations proportional to the rotational parameter a excepting the gravito- 
magnetic term - the second term in square brackets on the right hand side of ( 1B4I) can be 
neglected. These terms lead to terms in our final equation ( l48l) proportional to a product 
of two small parameters - oc a(/i/r)^. Thus, we can consider an effectively Schwarzschild 
problem with inclusion of only one term determined by rotation of the black hole. Addition- 
ally, the time derivatives of pi and can also be neglected since, similar to the Newtonian 
problem, they result in terms proportional to {td/ttw){h/r)'^ in the same equation. On the 



other hand, the time derivatives of the per turbed "horizonta. 



retained in a certain combination, see e.g. 



Papaloizou fc Lin 



" velocities and f should be 



(119951 ) ■ DI and below. The rea- 



son for this is the well known degeneracy of the Keplerian gravitational potential resulting 
in closed orbits of particles around a Newtonian source of gravity. Owing to this degeneracy 
a certain combination of terms proportional to the time derivatives (we call it later the 
"resonance combination") plays an important role at large distances from the black hole, 
where the gravitational potential is approximately the Newtonian one provided that a <^ 1. 



^ This is analogous to the freedom of choice of a gauge in Quantum Field Theory and General Relativity. 
^ Note that we as sume in this paper that ther e is no surface forces acting on the disc, for a more general case this requirement 

does not hold, see 



Ivanov fc Papaloizoul 1 2008l 'l for details. 
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3.2 The subset of equations describing the disc twist 

Taking into account the simpUfications discussed above let us consider equations of Appendix 

B. At first let us deal with the 'horizontal' part of the problem. It turns out that it is 

more instructive to use certain linear combinations of equations ( IBll) and ( ]B3[) instead of 

themselves. Neglecting the terms proportional to a and pi in these equations and taking the 

linear combination (jBl]) x j _ ^ x (2f/^) we remove the time derivative of v"^ 

from the result, thus obtaining 

1 IF d f mf\ TfP 

U'^P^ + JjjVfP^^^ + Kld^ y^'^TFJ = ^^"P^-^ + Wf^^^^"' - rWd.T:-). (45) 



In the Newtonian limit, when r — )■ cxd, equation (H5j) reduces to the continuity equation. 
Another equation reducing in the Newtonian limit to the ip component of perturbed Navier- 
Stokes equations can be obtained from (IBip and (1B3P taking another linear combination 
r( lB3l) X f/^ - (EI]) X t/^)/p with the resuh 

+ -^d^v^ + + V^' + 4r(9(rf - rWd^T^) = 0. (46) 

With our simplification being adopted equation (lB2p has the form 

^U^i,r + E^Q r _ 2^v^ + U(r:^ = Wr^. (47) 
Ki r ^ KiU'P p p 

It is easy to show that it reduces to the r component of perturbed Navier-Stokes equations 
when r — )■ oo. 

The vertical equation (lB4p plays a special role in our analysis. In fact, as we see later 
it most important for determination of the dynamics and form of twisted configurations of 
accretion discs. In the approximation we use, i.e. setting all terms proportional to a except- 
ing the gravitomagnetic term in ( IB4p . and, accordingly, in ( IB5tlB7P to zero we can bring 
this equation to the so-called divergent form, which refiects symmetries of the Schwarzschild 
space-time with respect to spacial rotations. This is analogous to the Newtonian problem, 
see PP and DI. For that we would like to express the density perturbation pi and the veloc- 
ity perturbation v'^ entering flB4p through v^' using equations (14511470 . Since the resonance 
combination of terms containing the time derivatives does not enter in the vertical equations 
we can set ■y'^ = v'" = when operating with equations (14511470 . After this being done we 
express the term containing d^pW in (H5l) through the velocity perturbations using equations 
(j27j) and ( 1471) . and substitute the result in ( H5i) . We express the velocity perturbation v'^ 



entering explicitly in fl45|) and flB4p . in terms of with help of fl46|) . and substitute it in 
these equations. Then, we use the relation obtained from ( H5!) to express the density per- 
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turbation through and the viscous terms and substitute the result in flB4l) . The resulting 
expression is integrated over ^, where we note that the term containing the derivative of the 
perturbed pressure, pi, vanishes after integration. This expression also contains derivatives 
over ^ in the viscous terms (]B5IIB6|) under the integrals. These are integrated by parts taking 
into account that the corresponding surface terms are equal to zero. After a rather tedious 
but straightforward calculation along the lines outlined above we obtain from ( ]B4p a rather 
simple equation 

L A2 J K2 

- J di { dri^rK.Kluya^v^ + T^K.KlT^^), (48) 

where we remind that S = f pd^ is the surface density and the bar stands for the quantities 
integrated over C,- 

In Appendix C we show that when a = equation (l48l) may be brought in a special 
divergent form determined by rotational symmetries of the Schwarzschild space-time. 

Provided a background model of a flat relativistic disc is specified equations ( H6M81) form 
a complete set. 



3.3 A model case of an accretion disc having isothermal density distribution 
in the vertical direction 

In order to transform the set of equations (1461148 p to a simpler form we should specify a 
dependence of density p on ^. Since an accurate treatment of the vertical distribution of the 
density is not important for our purposes we would like to consider the simplest possible 
case of an isothermal disc 

P = Pcexp(- — ), (49) 

where the central density pc and the disc half-height h are functions of r. 

Similar to the Newtonian problem it is easy to see that distributions of the velocity 
perturbations of the form 

v''^ = ^{Ai sirup + A2 cos (f) v'^ = $,{Bi sirup + B2 cos ip) (50) 

satisfy equations ( H6|l and ( H71) provided that the kinematic viscosity u does not alter with 
height and the amplitudes Ai, A2, Bi and B2 are assumed to be functions of the time and 
radial coordinates. 
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Introducing the complex notation 

A = A2 + iAi, B = B2 + iBi and W = + = /9e 
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(51) 



and using equations (|27|) , l^^ij , , (1B5P and (1B7P we can rewrite (H^ and (HTIl in the form 

-2 3 



respectively, where 



B 



{i - a)nB - 2nA = -{i + a)U'^nW', 



R- 



6 
R 



R-3 



A-(i- a)nA + ^B = -^iaK, (U^YV^WN', (52) 

(53) 

^ = R^iR-lyi- 

Note that k is the relativistic epicycl ic frequency for a slight l y perturbed circular orbit in 
the Schwarzschild space-time, see e.g. lAliev fc Galtsovl (Il98ll ): iKatd (Il990l ). 

As we have mentioned above equations ( l52l) . ( I53l) contain a degeneracy in the Newtonian 
limit r — 7- 00, for a formally inviscid disc, due to the well known Keplerian resonance between 
the mean and epicyclic motions. It is due to this resonance the " small" terms A and B should 
be retained in these equations. In order to clarify the physical meaning of this resonance let 
us consider the limit r — )• 00 and set a = in equations (152|) and (153|) . We get 



A-inA 



-B = 



(55) 



from equation (1521) . and 

B - iilB - 2fiA = -iW^n W, (56) 

where is assumed to have its Keplerian value, = r~^/^. Multiplying fl55|) by 2i and 
adding the result to (l56l) . we see that the terms proportional to Q on the left hand sides of 
( l55|) and ( l56|) cancel each other, and we obtain 

B + 2iA = -iU'^VL W (57) 

On the other hand we must have (|B|, |A|) ~ (trf/tf^)(|f2A|, |f^B|), where we remind that 
td ~ is the characteristic dynamical time scale, ttw is the characteristic time scale of 
evolution of twisted disturbances and t^/tty^ ^ 1. Taking this fact into account we see that 
the term on the right hand side of fl57|) is of order of {tii/ttw)\^A\. Therefore, it can be 
neglected in the leading approximation in (l56ll . and we get either from this equation or from 
( l55l) an algebraic relation between A and B: 

B = -2iA, (58) 

which is valid only in the leading order in ta/tt^. Thus, in this order the perturbed velocities 
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can be described by only one complex amplitude, either A or Oo- Let us show that the same 
property approximately holds for the general relativistic problem as well. For that let us 
consider a linear combination of (l52l) and (l53l) with coefficients chosen in such a way that 
the term proportional to A is eliminated 



B - 



2Q 



[i — a)VL 



A 



1 + 



K 



{i - a)2fi2 



3ia 



I — a 



W. 

(59) 



The combination on the left hand side plays a role in dynamics of the disc only when we 
consider a non- relativistic, low viscosity disc with a < h/r while in the opposite limits it 
can be set to zero. Therefore, we can consider all background quantities in this combination 
equal to their Newtonian limits, set a = there and use equation fl58l) to express A through 
B on the right hand side. We get 



B 



1 + 



{i - a)nB 



(z + a)f/'^n- 



3ia 



a 



(60) 



Equation ( l60l) describes the horizontal part of the problem in the leading order in td/Uwi for 
all possible values of a and r. It is one of two final equations describing a fully relativistic 
twisted disc around a slowly rotating black hole. 

In order to get the second final equation describing our problem let us consider the 
vertical equation ( l48l) . At first we use equations ( HOllSTl) and (IB 71) to rewrite it in terms of 
the complex variables and perform integration over ^. After this being done we get expression 
proportional to fj and S/i^ under the derivative over r. We substitute f p3|) and fl35|) and take 
the constant factors out of the derivative. Then, we use equation (!28l) to express U'^ in terms 
of M and the explicit form of T^'^ given by fl30|) . After the resulting expression is divided over 
S one can see that it contains M and S only in the combination M jH. This combination is 
expressed in terms oi 5 = h/r with help of ( l36l) . We finally obtain 



2 A2 V ^ J 

2^1/2^^3/2^ <9r \ ^ KlU^L 
Equations (j60ll and (j6Tl) form a complete set. 



(i + a)B + af/'^W ) y . (61) 



* As discussed in e.g. DI the relation l|58|l follows from the fact that slightly perturbed circular orbits of free particles in the 
Newtonian potential are ellipses with a small eccentricity. 
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4 A STATIONARY TWISTED DISC 



Non-stationary solutions of equations ( 1601) and (1611) will be discussed in a separate publica- 
tion. In this paper we would like to consider a stationary problem setting time derivatives in 
(|60|) and (|6T]) to zero. Then, we express B through W using (|60|) and substitute the result 
in (16T]) to obtain 

where * stands for the complex conjugate, 



/(a, R) = {l + a'- ?,raKl) ^ + a, (63) 

ait(a -I- 2z) — 6 

we use the coordinate R as an independent variable, and the explicit expression for the 
gravitomagnetic frequency VL^t (IBlip is used. It is easy to see that in the non-relativistic 
limit -R — )■ oo and i^i — )■ 1 the function f{a , R) tends to the function /(a) - 



{a+2i) 

introduced in paper iKumar fc Pringld ( ll985l ). hereafter KP. 

A character of solutions to equation depends on only two independent parameters - a 
and 

S = 6J^/m (64) 



4.1 Limiting cases 

A number of authors considered equations describing a shape of a stationary twisted disc far 
from a rotating black hole. They used either a purely Newtonian analysis with the gravito- 
magnetic term treated as originating from a Newtonian force responsible for precession of the 
disc rings (e.g. PP, KP) or took into account certain post Newtonian tern is, which may play a 



Ivanov &: lUarionov 



significant role in de termining of the stationary disc c onfigurations (e.g. 
(Il997l ). hereafter II; iLubow. Ogilvie fc Pringld (120021 )). Obviously, equation ( l62l) should be 



reduced to equations exploited by previous authors in the respective limits. 

In particular, in the paper II a low viscosity (a ^ 1) stationary twisted disc was con- 
sidered and a simple equation describing the stationary configurations was derived taking 
into account only the most important post-Newtonian correction to the first term on the 
right hand side of ( l62l) . In the approximation used by II we neglect the second term on the 
right hand side of ( l62l) . which is proportional to a and set the values of all quantities there 
equal to their Newtonian limits (i.e. = 1, If^ = D = 1 and Ki = 1) and assume 

that a = everywhere excepting the so-called 'resonance' denominator in the first term in 
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(!63|) . which should contain a term proportional to a and a first relativistic correction. In 
this approximation we have /(a, R) = ^jrrKT obtain from equation ( l62i) 

d 1 fw+^^W = 0, (65) 



c/xi (1 — Sia^^xf) dxi S'^ 
where Xi = This equation coincides with equation (33) of II. 

PP derived dynamical equations describing a Newtonian twisted accretion disc in the 
limit of relatively large viscosity a ^ 6. Stationary solutions to these equations were ex- 
amined in the paper KP, for a disc with isothermal distribution of density with height. In 
order to obtain a fully Newtonian limit of fl62|) and compare it to equation used by KP we 
set If^ = 1, If^ = _R~^/^, and i^i = 1 in (162!) and (!63|) . We also introduce the independent 



variable X2 = y/R^xi = a/ Rms/R, use the form oi D = 1 — X2 employed by KP and 
parametrisation of 5, 5 = {H^/ Rms)x2^^^ ^^Vd, where g is a parameter as in KP. In this 
way we obtain from ( 162|1 equation studied by KP 



xr^^(r(«)A 

(1X2 



X2)^^ 
dX2 



6a Aw ) + 16mx^^^W = 0. (66) 
dxo / Hi 



4.2 An almost inviscid twisted disc 



An important limiting case of (1621) is obtained by formally setting a = in this equation. 
When the disc gas rotates in the same sense as the black hole and, accordingly, a > the 
behaviour of the disc inclination angle differs drastically from the standard picture of the 
well known Bardeen-Petterson effect, where it is expected that the disc aligns with the black 
hole equatorial plane at sufficiently small radii. Contrary to this picture when considering 
the shape of low viscosity twisted disc at large distances R ^ R^s II found a phenomenon 
of spatial oscillations of the disc inclination angle /3 for the case of a > while in the 
opposite case the disc does align with the black hole equatorial plane. Our fully relativistic 
analysis allows us to consider this effect in more details and show that the amplitude of the 
oscillations can be significantly amplified close to the last stable orbit. 
Setting a = in (1^ we get 

where 

_ R^/^D _ 2AaD 
^-1^' ^-^KfU^^- ^^^^ 

Note that equation (167|) depends parametrically only on 6. Since this equation has only real 
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coefficients its solution may be taken to be real. Thus, in this limit only the angle /3 changes 
with radius while 7 stays unchanged. 

It is also important to note that the influence of the black hole gravitomagnetic force on 
the disc shape is most significant when 6 <^ 1. Also, it is most natural to suppose that the 
rotation parameter a > 0. The case of 5 ^ 1 and a > is considered in more details below. 
Note, however, that the bulk of our expressions below are formally valid for the case a < 
as well. Although in the case considered below we may choose W = /3 we continue to use 
W to keep our notation uniform. 

4-2.1 The shape of an almost inviscid stationary twisted disc close to Rms 

At ffist let us discuss the form of solutions to (1671) close to the last stable orbit assuming 
that X = R — Rms ^ 1- For that we change the independent variable to x in ( l67|) and 
assume that all quantities entering this equation, which have finite values at Rms are equal 
to these values. Thus, we take into account only the dependencies of D and 5 on x in such 
approximation. We use equation fHOl) . which tells that D ^ x^ jll close to Rms- 

We parametrise the dependency of 5 on x as 5 = bms^'^'' ■, with bms being proportional 
to We are going to consider only the values of e < 1/2, which physically corresponds to 
the condition that the surface density S tends formally to zero when R — )■ Rms- Equations 
( 1371) and ( l38l) tell that e = 1/5, bms = 0.375=,, and e = 3/20, 6ms = 0.515=,,, for the opacity 
dominated by the Thomson scattering and the free-free processes, respectively. 

In this way we get from ( 167|) 

i_a;2 _1 w + xx~^'y^ = 0, (69) 
x'^ dx dx 



where 



^ 24a[/" 



'ms 



is evaluated at Rms- Using the parameter 5 given by flMj) we can rewrite f l70|) in the form 
X = X^'"^) where xtk = 8.4 ■ 10~^ and Xff = 4.4 ■ 10~^, the indices Th and // stand for 
the cases of the Thomson and free-free opacities, respectively. When a > the solution to 
equation regular at R = Rms may be expressed through the Bessel function, Juiz), as 

W = Cx-i/Vi/2(i-2.)(^), (71) 

where 

-x'-^\ (72) 



;i-2e)- 
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It is easy to see that the solution (I7T]) tends to a non-zero constant when z ^ 0. This fact 
may be used to express the constant C through the value of W at x = 0, Wq, as 

r(Jp^)(^y W,, (73) 



,2(1 - 2e) J \2{1 - 2e) 

where T{x) is the gamma function and we use the well known approximate expression for 
Ju{z) at small z. Since the solution contains the Bessel functions of a real argument it 
describes oscillations of the inclination angle /£. 



Let us estimate a change of the amplitude of these oscillations over the range < x < 1, 
in the limit of 5 ^ 1. As follows from flTOj) and f l72|) in this limit the variable z takes values 
^ 1 even for x < 1 and, therefore, the region can be divided on two domains: z < 1 and 
z > 1. In the domain z < 1 the value of W monotonically decreases with x, and, typically, 
it decreases in six times before oscillations of the inclination angle begin. Thus, we have 
Wo/W(z ~ 1) ~ 6. When z ^ 1 we use the asymptotic expansion of the Bessel function in 
powers of z~^ simplifying f lTT]) to 



W^CJ^cosiz- -^-1-) . (74) 

\ TTxz \ 21 - 2eJ ^ ^ 

Assuming that this expression is approximately valid even when x ~ 1 or 2; ~ 1 and that 

I — 2e ~ 0(1), we can roughly estimate 

W(X ~ 1) ~ X'(l-^)/2(l-2e)^(^ _ ^75^ 

Using the quantities S and x iii dZH]) and remembering that Wo/W(z ~ 1) ~ 6 we obtain 

Wo/W(x ~ 1) ~ ( Q^(^-W-2e) ) ^-(l-.)/(l-2.)_ (7g) 

Since Gx*-^"^-''^^''^"^'^'' ~ 0(1) for the parameters we use we set this factor to unity in our 
expressions below. 

From equation (1761) it follows that when a — )■ and 5 ^ 1 a rise of the amplitude of the 
oscillation with decrease of x near the marginally stable orbit can be quite large. 

4-2.2 The shape of the disc at large radii: the case of a > and 5^1 

II showed that when the limit 5 ^ 1 is considered a low viscosity disc inclined at large 
distances with respect to the equatorial plane gets twisted at a characteristic scale 

R2 ~ 5-'/' > Rms. (77) 



^ Let us stress again that in the case of a < the solution contains the Bessel function of an imaginary argument, which 
describes a monotonic decrease of the inclination angle of the disc with decrease of x. 
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At such large scales we can use equation (!65|) considering the limit a — )■ there. In this limit 
we get 

d^ d ~ 

a;i-^W-2— W + 96r\fW = 0, (78) 
dxf dxi 

where we remind that xi = and that we consider the case a > 0. Note that following 

II we neglect a very slow dependence of 5 on i? for R ^ 1 and set 5 = 5* in (ITHj) . The 
solution to ( ITHI) can be expressed through the Bessel functions 

W = xf\A,J_,/,{z,) + A2J3/5(^l)), (79) 

where 

= -V6S-'xl^\ (80) 
l— I 5 

and Ai and A2 are constants lj. The limit of large R^ R2 corresponds to the limit of small 
Zi. In this limit the first and second terms in the brackets in fl79|l being multiplied by Xi^'^ 
tend to a non-zero constant and to zero, respectively. This fact may be used to express the 

constant Ai through the asymptotic value of W at infinity, Wgo, as 

/ 5 \3/5 ^3/5 

In the opposite limit oi R <^ R2 we get a simple approximate expression for W(i?) from 
(179]) and (ISTi) using the asymptotic expressions for the Bessel functions at large values of 
their arguments 



W^^l^^R-y 



cos (2:1 + ^) + sin [z,~^y (82) 

Equation f lS^ tells that in the region R^g <^ R <^ R2 the inclination angle oscillates with 
the amplitude of oscillations proportional to R^^^^ . 



4.2.3 A WKBJ analysis of (6^ in the limit o/(5 < 1 

The asymptotic expression fl7^ and (!82|) can be matched together by a WKBJ solution. 
Indeed, in the limit 5^1 the ratio of the quantities A and b entering (|67|) . A = A/6, is 
large for the range of R such as z and Zi defined in equations (172]) and (17^ , respectively, are 
large. Applying the standard WKBJ scheme we find the solution in this region in the form 

W ^ cos (^j^ VldR + <i)wKB^ , (83) 



Note that an analogous equation of II contains a misprint. Namely, the power of their variable y2 (which is equal to our 
variable zi) in front of the square brackets in their equation (37) should be 3/5 and not 3/2 as in the text. 
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where the constants Ci and (pwKBj must be chosen in such a way that the solution (!7il) is 
matched in the corresponding asymptotic limit. It is easy to see that in order to have such 
a matching we should require that 



TT 1 — e 



mKBj 



2 1 - 2e' 
and 



(84) 



where it is assumed that Ki and U'^ are evaluated at = Rms and we use the facts that 
Rms = 6 and that D ^ 112 close to Rms- 

In the limit i? — ?■ oo we can set the values of A and h in front of the cosine in ( !83l) equal 
to their Newtonian values. The integral in ( 183|) can be represented as l[R) = \/\dR = 
I — ^/l\dR, where / = \/xdR. Taking into account that the Newtonian value of 
A = 245-2/2-9/4 we have VxdR ^ ^6-^R-^/^, and, accordingly, 

W ^ cos (^^~6-'R-'/' - I - c^wKB^j . (86) 

The expression fl5B|) can be matched to the asymptotic expression (1521) provided that the 
constant Ai and A2 are appropriately chosen. A simple calculation gives 



Ai = )J^Ci cos (^/ + 



W™-^)/cos^, 



A2 = sin (/ + (pwKBj + ^) / (87) 

Equations fl7ip . f l79p and 083 p provide expressions for an approximate shape of a low 
viscosity twisted disc in the whole allowed range of R, Rms < R < 00. Note that although 
we derive these equations assuming that a > they are formally valid for negative values of 
a as well. 

Equations fl73|) . flHTl) . fl85|l and flHTj) allow us to relate the asymptotic value at large radii. 
Woo, to the value of W at the last stable orbit, Wq, as 

Woo = Ctot{S)Wo, (88) 
where an explicit form for Ctot{S) follows from these equations. It is important to note that 
Ctoti^) oc cos(/ + (pwKBj ^ ^) as follows from equations flHTl) and (1571) . Thus, for values of 
5 such that cos(/ + 4>wkbj ^ i^) = 0; = while Wq 7^ 0. This describes a peculiar 
"resonant" solution to equation ( 1671) . where a regular at the last stable orbit solution is 
matched precisely to the solution of fl78|) proportional to A2 while Ai = 0, see equation fl79|) . 
From equation (!68|) it follows that the integral / in the cosine can be written in the form 
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I = 6^^ I, where / does not depend on 6, we have Ixh = 1-01 and /// = 0.69. This allows us 
to write the condition for obtaining the resonant solutions as a condition that the parameter 
6 has discrete values 

where k is an integer. Thus, a stationary twisted disc described by such a solution has its in- 
clination angle going to zero with radius. Note, however, that inclusion of effects determined 
by the viscosity coefficient a leads to a partial suppression of this resonance. As we see later 
numerical results show that when the resonant condition (l89i) is fulfilled and a small but 
nonzero value of a is considered the value of Ctot at the resonant values of 5 is much smaller 
than its neighbouring values instead of being equal precisely to zero. 

Neglecting the possibility of the resonances we can roughly estimate the dependence 
of Ctot on 5 using very simple arguments. Indeed, from equation ( l82l) it follows that the 
inclination angle scales with as W oc R^^^^ for radii Rms < R < R2, where R2 is given 
by equation ( I77|) while when R> R2 it stays approximately constant. Thus, we can roughly 
estimate W(i? ~ Rms) ~ R2^^^oo = ^~^^^^^oo- This estimate does not take into account 
the growth of the amplitude of the oscillations close to the last stable orbit, which is described 
by equation ( 1751) . Using this equation and assuming, for simplicity, that the numerical factor 
g^(i-e)/2(i-2e) _ ^ ^^leYe, wc estimate Wo/W(i? ~ Rms) ~ Combining both 

estimates we get 

Ctot ~ 5(i-)/(i-2^)+Vio. (90) 



4.2.4 Restrictions 

We derive our equations assuming that the inclination angle (3 is small. Also, we suppose 
that the NT model of the flat background disc with a constant value of the parameter a is 
valid. These assumptions put rather severe limitations on validity of our analysis, especially 
in the case of the low viscosity twisted disc considered in Section 4.2. 

Since in this case 7^0 and /3 ~ W the condition /3 < 1 everywhere in the region 
R > Rms leads to the condition Wq < 1. Equation (pO|) tells that in this case the asymptotic 
value of the inclination angle at spacial infinity should be quite small 

< W,„, = 5(l-)/(l-2e)+l/10. (91) 

Another somewhat more stringent constraint stems from the fact that the shear veloc- 
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ities induced in the disc by the disc twist should not be too large. When considering the 
oscillations of the disc inclination angle at scales ~ R2 II supposed that a characteristic 
amplitude of the shear velocities should not exceed the sound speed, Cg. In the opposite 
case the shearing instability could make the disc more turbulent, thus increa sing an effective 



value of a. Another possibility is the presence of shocks in the disc (e.g. [Fragile &: Blaes 



(120081 )). which may heat up the disc thus effectively increasing the relative disc thickness 5. 
II obtained a condition that the relativistic oscillations do not induce such velocity 

4/5 

amplitudes at the scale ~ R2 provided that Woo < 25* . Note, however, that our numerical 
results suggest that viscosity effects suppress the oscillation amplitude, and, accordingly, 
the velocity amplitude more efficiently than it was assumed in a simplified analysis of II. 
Therefore, the constraint of II may to be too restrictive for more realistic models taking into 
account a small but non-zero value of a, see the next Section. 

In order to estimate a characteristic amplitude of the shear velocities induced in the disc 
we assume that it is of the order of Vsh = h'B, where B is defined in equation ( 1511) and we 
remind that it can be considered as a real quantity as far as equation (|67j) is concerned. 
From equation (l60l) it follows that B ^ [/"^W when B and a are equal to zero there and 
K ^ 0, which is appropriate for the radii close to Rms- On the other hand, from the results 
of Section 2.4 it follows that = Lf^h/R. Introducing the variable v = Vsh/cs we have 
V = RW ^ 6W'. From equation ( 1741) it follows that when x < 1 but 2; > 1 we have 
W ~ ~ and from equation ([75]) it follows that C ~ X~(^~^')/^Wo, 

and, therefore, 

V ~ 6Wox'/^a;-(^+'\ (92) 

Since x = X^~'^ and the factor x'^^'^ ~ 1 we can roughly estimate that x ~ ^ "^ ( 1921) . We 
see that the condition {; < 1 is fulfilled when 

X > x,rit = (6Wo)^/(^+^)r^/(^+^). (93) 
According to our criterion our analysis is valid close to Rms provided that we demand that 
Xcrit < 1, at least. This leads to inequality Wq < 5^/6, which can be reformulated in terms 
of Woo using equation ( pOj) 

Woo < W,„i576, (94) 

where Wcrit is defined in ( 19T1) . We see that this inequality is stronger than the inequality 
( 19T|) although the difference is not quite significant since the parameter e is typically rather 
small. 
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Figure 1. Left panel: The dependence of /3 on the radial coordinate R is shown. Since the problem is a linear one we set 
hereafter the value of the inclination angle at large radii, /3oo, equal to unity in all Figs. The parameter S = 10~^. The solid 
curve is obtained by setting o = in I I62I I. The short dashed curve represents our analytical results discussed above. We 
use the expression for /? given by equation (I83I I in an intermediate region of R, when R — Rms <^ 1 the expression 1711 1 is 
used and when R ~ R2 or larger we use equation I l79t . The dot-dashed, dot-double dashed, long dashed curves correspond to 
a = 10""*, 10~^ and 10"'^, respectively. Right panel: The result of integration of equation II62I I presented in parametric form. 

= i?eW = /3cos7 and ^2 = ImW = /3sin7, 5 = 10~^ and a = 10"'^ for this plot. The circle, square, triangle (pointing 
up), triangle (pointing down) and the star correspond to i? = 6, 6.1, 6.2, 6.3 and 7, respectively. 

4.3 Numerical results 

We integrate numerically equation ( l62ll starting from the last stable orbit. Similar to the 
case of almost inviscid twisted disc considered above there is a partial solution regular at 
the last stable orbit and we use this solution to specify the boundary condition at Rms- 

We consider below the case of Thomson opacity only. The case of free-free opacity is 
quite similar. 



4-3.1 Dependencies of P and 'y on R 

At first let us discuss the dependency of /3 on i? for the case of relatively small viscosity 
parameter a ^ 10^^ and a > 0. This case is represented in Figs. [1] (left panel) and [21 for 
6 = 10~^ and 0.1, respectively, lj As seen from these Figs, the oscillations of the inclination 
angle are quite prominent for the case when a = in equation ( l62|) as well as for the case 
of quite small a = 10~^. Also, it is seen that our analytic theory described above gives the 
curve (the left panel in Fig. [T]), which is in a quite good agreement with the numerical one 
calculated for a = 0. Effects of viscosity tend to smooth out the oscillations, which disappear 
for the curves with a = 10^^ and 10^^ on the left panel in Fig. [H The similar tendency is 



Note that we show the ratios /3(i?)//3oo in these Figs, using the fact that our problem is a linear one. This ratio can be 
arbitrary large provided that /3oo is sufficiently small. 
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Figure 2. The same as on the left panel in Fig. [T] but for 5 = 10^^. The solid, dotted, short dashed, dot-dashed, dot-double 
dashed and long dashed curves correspond to a = 0, 10~*, 10^'^, 10~^, 10"^ and 1, respectively. 

seen in Fig. [2] where, additionally, the case of relatively large a > 10^^ is shown as well. 
It is clear from these Figs, that when a is small the inclination angle grows with decrease 
of R either on average or monotonically. Thus, in these cases the disc does not align with 
the black hole equatorial plane and the Bardeen-Petterson effect is absent. Note that the 
smooth curves do not have sharp gradients of the inclination angle and, accordingly, they 
correspond to solutions, where the shear velocities induced in the disc are relatively small. 
Therefore, these solutions may be physically realistic even for relatively large values of the 
inclination angle at large radii, (3oo- 

Instead of showing dependencies of the second Euler angle 7 on i? we plot curves corre- 
sponding to our solutions on the plane (\E'i = /3cos7, \E'2 = /3sin7), which depend paramet- 
rically on R. Since the angle 7 does not change significantly with R for very small values 
of a and the corresponding curves are close to straight lines, the case of small viscosity 
parameters is represented by only one curve corresponding to the solution with a = 10~^, 
see the right panel in Fig. [H One can see from this Fig. that this curve has a form of a tight 
spiral. In this case /3o ~ 8 and 7 = 0, so the curve starts at i? = R^s when ~ 8 and 
\E'2 = 0. It spirals clockwise with increase of R, i.e. in the direction opposite to the orbital 
motion and black hole rotation, which are counterclockwise. 

The dependence of /3 on i? as well as the solutions in the parametric form for 6 = 10^^ 
and larger values of a = 0.05, 0.1, 0.2 and 1 are shown in Fig. [31 From the left panel in Fig. 
|3] one can see that in this case the inclination angle monotonically decrease in practically 
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Figure 3. Same as in Fig.[T]but for larger values of a. The solid, dashed, dot-dashed and dot-double dashed curves are plotted 
for a = 0.05, 0.1, 0.2 and 1, respectively. Curves of the same style correspond to the same a on both panels. Different symbols 
on the right panel correspond to different values of R, R = 100, lO'' and 10'' for the circles, squares and triangles, respectively. 



the whole range of R while there is a small increase of the angle close to Rms- In general, the 
ratio of Po/Poo is quite small and this may be interpreted as manifestation of the Bardeen- 
Petterson effect. Contrary to the case of the low viscosity discussed above the initial value 
of is smaller than one, thus for the curves on the right panel of Fig. [3] /3 = a/^&i + 
increases its value when the curves spiral from R^s to larger R. The direction of spiralling 
with the increase of R is also clockwise similar to what is shown on right panel of Fig. [1] 

However, returning to Fig. [2] we see that for 6 = 10~^ the Bardeen-Petterson effect is 
absent even when a ~ 1 in the sense that the ratio (3q / /3oo remains to be of order of unity. 

Thus, in the case of a > the Bardeen-Petterson effect occurs only when the viscosity 
parameter a is sufficiently large and 6 is sufficiently small. 

The case of a < is shown in Fig. m for 6 = 10^^ and a range of values of a. As is seen 
from this Fig. contrary to the previous case the disc always aligns with the black hole. This 
is valid even in the case of setting a = in fl62|) . which is described by the solid curve in 
Fig. m The curves with small values of a are close to this limiting curve. It is interesting 
to point out that as seen from the right panel in Fig. H] the direction of spiralling when R 
runs from smaller to larger radii is opposite to the case of a > 0. Now the curves spirals in 
the direction of the orbital motion but, again, opposite to the direction of the black hole 
rotation. Thus, the direction of spiralling is always opposite to the black hole rotation. This 
may have some observational consequences. 
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Figure 4. The same as in Fig. [T] but for tlie case of tlie black liole rotating in the opposite sense with respect to the disc 
rotation, a < 0. The solid, dotted, short-dashed, dot-dashed, dot-double dashed and the long-dashed curves correspond to 
a = 0, 0.01, 0.05, 0.1, 0.2 and 1, respectively. Curves of the same style correspond to the same a on both panels and positions 
of different symbols on the right panel correspond to the same R as in Fig. \3\ 

4-3.2 An analysis of behaviour of the solutions for the case a > on the parameter plane 
(a, 6) 

As we discuss above our equation contains two independent parameters - a and S. It is of 
interest to describe how the most important characteristics of the solutions such as the ratio 
Pol Poo depend on them. 

In Fig. [5] we show the dependence of this ratio on 5 for various sufficiently small values 
of a. One can see that the curve corresponding to a = experiences a quite non-monotonic 
behaviour with series of peaks at some successive values of 5. Since the analytical curve 
based on our WKBJ theory developed above is in excellent agreement with the numerical 
one we can identify these peaks with the resonant solutions of (!62|) with a = for which 
the value of the inclination angle tends to zero when i? — )■ oo. The positions of these peaks 
are, accordingly, given by equation (!89l) . Note also that our simple estimate (!90!) of the 
dependence of Po/Poo on 6 provides a very good approximation to the numerical curve when 
the resonance peaks are not taken into account. 

From the analytical theory it follows that when a = formally the amplitude of peaks 
in ([5]) is infinite. When a small value of a is taken into account the peak's amplitudes get 
finite values decreasing with increase of a and/or 6. For a relatively large value of a = 10~^ 
the peaks practically disappear and the corresponding curve strongly deviates from the 
theoretical one. In the case the values of Po/Poo ~ 10 in the whole range of considered values 
of 6. 
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Figure 5. The dependence of the ratio of the inchnation angle at the last stable orbit to the inclination angle at a large radius, 
/3o//3oo, on the parameter S. The solid curve is calculated for a set to zero in II62I I. the short dashed curve represents C^J^{S), 
where Ctot is defined in equation l|88|l while the dotted curve represents our simple estimate l|90|l . which does not take into 
account the effect of the resonances. The dot-dashed, dot-double dashed and the long dashed curves are for a = 10~*, 10~^ 
and 10"'^, respectively. 
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Figure 6. Levels of constant ratios of /3o//3oo on the plane {5, a) of the parameters of the problem. Arrows show the values of 
/9o//3oo for each curve. The dashed curve on the right hand side of the plot separates a region on the parameter plane, where a 
change of /3 with R is larger than 10 per cent of /3oo (to the left of this curve) from the region, where the black hole influence 
on the disc twist is relatively small and the change is less than 10 per cent (to the right of the curve). 



In Fig. |6]we show the curves of constant values of the ratio of Po/Poo on the plane {5, 
a) of the parameters entering equation fl62|) . for the case of a > and the value of 6 in the 
range 10"'^ < 5 < 10. When this ratio is small the disc aligns with the black hole equatorial 
plane at small radii. One can see that the area of small values of /3o//3oo is situated at the 
upper left part of the plot. The opposite case, where (3^/(3^ > 1 corresponds to a situation 
when the Bardeen-Petterson effect is absent. For the considered range of S the area in the 
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Figure 7. Dependencies of /9 on along tlie curve in Fig. [6] wliere fio/ Poo = 1- The values of 5 are equal to 10 "^,10 ^,10 ^ 
and 1 for the solid, dashed, dot-dashed and dot-double dashed curves, respectively. 



plot, where Po/Poo > 1 is larger than the area corresponding to the small values of this ratio. 
The curves corresponding to the levels of /3o//3oo ^ 10 exhibit oscillatory behaviour. This 
effect is most probably determined by the resonant behaviour of Pq / Poo at small values of a 
considered above. 

The dashed curve on the right hand side of the plot separates the regions on the parameter 
plane where the quantity \P{R)/Poo — 1| is smaller /larger than 0.1, for all values of R. The 
latter condition is satisfied to the left of this curve, where the disc twist induced by the black 
hole is sufficiently large. Thus, the black hole significantly infiuences on the disc only when 
the value of 6 is sufficiently small, 5 < 0.1 — 1 depending on value of a. 

There is a curve in the plot corresponding to the level of Po/Poo = 1- The dependencies 
of P{R) with parameters a and 6 taken along this curve are shown in Fig. [71 As is seen 
from this Fig. when S is sufficiently small the angle P behaves in a non-monotonic way, 
though the curves are smooth and no oscillations of the inclination angle are observed. Such 
a behaviour looks realistic and it may have serious consequences for the disc structure itself 
since the inner parts of the disc strongly irradiate the outer part in such cases. Also, it may 
have some interesting outcomes for spectra modelling, etc. 



5 DISCUSSION 



In this Paper we derive dynamical equations describing evolution and stationary configura- 
tions of a fully relativistic thin twisted disc around a slowly rotating black hole. We assume 
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that the disc inchnation angle /3 ^ 1, the black hole rotational parameter a is small and 
that the background flat accretion disc is described by the simplest Novikov-Thorne model 
with a constant value of the Shakura-Sunyaev parameter a. With these assumptions being 
adopted a characteristic evolution time scale of such a disc is much larger than the dynamical 
time scale. This helps to simplify our analysis to a great extent. It is shown that the final 
dynamical equations may be formulated as equations for two dynamical complex variables 
W describing the disc twist and warp and B, which describes shear velocities induced in 
such a disc. They may be brought to a form similar to what has been used for description 
of the classical Newtonian twisted discs, see equations ( 160|) and ( 16T|) . 

Since all our variables are well defined in the black hole space-time, the results of inte- 
gration of our dynamical equations may be directly applied to the problems connected with 
observations, e.g. obtaining time-dependant luminosity and spectrum of a twisted accretion 
disc, its shape as seen from large distances, etc.. 

We analyse stationary configurations of twisted discs in our model and show that they can 
be fully described by only two parameters - the viscosity parameter a and by the parameter 
6 = 5^/ \/\a\, where 5* is a characteristic opening angle of a background flat accretion disc. 

We consider in details the case of a very low viscosity twisted disc formally setting a = 
in equation f l62|) describing the stationary configurations and develop an analytic approach 
to the problem of finding of these configurations for the most interesting case 5^1 and 
a > 0, which is in excellent agreement with results of numerical integration of ( l62l) . We find 
that the oscillations of the disc inclination angle found by II in their model of a low viscosity 
twisted disc far from the black hole can be in a different regime in the relativistic region 
close to the last stable orbit. Namely, for the background fiat disc models having singular 
profiles with the surface density being formally zero at Rms as in the gas pressure dominated 
NT models the amplitude of disc oscillations and the radial frequency can grow indefinitely 
large when 5 — )■ 0. Additionally, in the formal limit a = there are "resonant" solutions 
occurring at some discrete values of 6, which have a finite value of the disc inclination angle 
at i? = Rms and the disc inclination tending to zero when i? — > oo. 

Effects of non-zero viscosity tend to smooth out the relativistic oscillations to an extent 
more pro minent than was obtained in an oversimplified model of II. This may explain the 



results of 



iNelson fc Papaloizod ( l2000l ) who didn't find prominent oscillations of the disc 



inclination angle in their SPH simulations of the Bardeen-Petterson effect. When viscosity 
parameter a gets sufficiently large (say, a ~ 10~^ — 10~^ for 6 = 10~^) the disc oscillations 
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are not observed. Instead of it there is a growth of the inchnation angle towards the last 
stable orbit, which is opposite to what is expected in the standard picture of the Bardeen- 
Petterson effect. When a gets even larger (say, 10~^ for 6 = 10~^) the inclination angle 
decreases towards Rms as in the standard picture. 

In general, there are three possible quahtatively different shapes of the twisted disc 
depending on values of the parameters of the problem. 

1) When a is very small the disc inclination angle oscillates, with the oscillation amplitude 
growing dramatically towards -Rms- When the disc inclination angle at large radii is not very 
small such a regime is probably unphysical. As we discuss in the text in such a case the 
shear velocities induced in the disc can easily exceed the speed of sound. This may either 
make the disc more turbulent thus effectively increasing the value of a, which may, in its 
turn, damp the oscillations. On the other hand this can lead to presence of shocks in the disc 
close to Rms- The energy dissipated in the shocks may make the disc thicker. This opens up 
a possibility of explaining of the presence of more compatible with observations thick discs 
close to Rms by such a mechanism. 

2) When a is moderately small the oscillations are absent. However, the Bardeen- 
Petterson effect is absent as well. In such a regime the disc inclination can either monoton- 
ically grow towards the last stable orbit or exhibit a non-monotonic behaviour decreasing 
with decrease of R at large radii and increasing at smaller radii, see Figs. [Hand [71 Since the 
corresponding curves are quite smooth we do not expect large shear velocities, which are 
proportional to the gradient of W. Therefore, such a regime may be physically allowed. This 
can lead to interesting possibilities of modifying the disc structure by strong irradiation of 
outer parts of the disc by inner parts, can be very important for modelling of spectra of the 
discs, etc.. 

3) Finally, when a is sufficiently large the Bardeen-Petterson effect is observed. The 
region of the parameter plane, where such effect is possible is shown in Fig. [6l for the 
considered range 10"'^ < S < 10. 

Although we consider the rotational parameter of the black hole as a small parameter of 
the problem we believe that the outlined picture remains qualitatively correct for moderate 
values of the rotational parameter as well, say when a ~ 0.5. 

We discuss possible generalisations of our results. At first one can solve time-dependant 
equations f l60|) and f[6T|) to clarify dynamics of the twisted discs in the relativistic regime. 
This will be done in a separate publication. Secondly, one can try to generalise the results to 
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the case of large values of the rotational parameter, a ~ 1. We expect that time-dependent 
equations will be rather cumbersome since in this case the evolution time scale will be of 
the order of the dynamical one and many terms neglected in our analysis must be taken into 
account. The generalisation to the stationary case looks relatively straightforward. However, 
it may happen that our twisted coordinate system based on the cylindrical coordinate system 
is less convenient when a ~ 1 than a twisted coordinate sy stem based on the spherical 



coordinates similar to what was introduced by lOgilvid (Il999[ ) for the classical Newtonian 
twisted discs. One may also try to include (at least, at some phenomeno logical level) the 
effects of back reaction on the background flat disc model determined by the large shear 
velocities induced in the twisted disc, and the disc self-irradiation. Additionally, as we discuss 
in Section 2.4 equations describing the flat disc models used in this paper are invalid very 
close to the last stable orbit, where the radial drift velocity gets larger than the sound 
speed, see equation (142|) . One may try to use a slim disc model close to Rms as a background 
model for our equations to overcome this difficulty. Finally, one may use the shape of the 
stationary configurations calculated in this paper to model different observed characteristics 
of particular sources. 
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APPENDIX A: THE ONE FORMS AND CONNECTION COEFFICIENTS 



Al The forms 



The transformation law may be considered as a rotation to another Cartesian coordinate 
system [t = t,xi = r cosip, yi = r sin tp, z\ = ^) with help of the rotational matrix entering 
([7]). Let us call this matrix as A. The one forms (11011131) are related to the forms associated 
with the new coordinates {T,xi,yi, zi) as the forms ([6]) are related to the forms associated 
with the old one {t,x,y,z). Explicitly, the transformation between the forms ([6]) and the 
forms associated with the coordinates {t, x, y, z) may be written as 



(Al) 



where uj are vectors having the corresponding one forms as their components, e.g. has 
the components etc.. B is a rotational matrix with components 



(A2) 



/ 1 \ 

cos(f) — sin 

sin(j) coscf) 

\ 1 / 

Analogously, the transformation between the vector containing the forms fll0m3p . uitw) and 
the vector containing the forms corresponding to the new Cartesian system, (jJcarNew^ may 
be written as 



1^ Car New 



Cw 



(A3) 



where the rotational matrix C is obtained from B by change of the angle: (/) ^ i/j. Accord- 
ingly, we have 



(A4) 



(A5) 



The rotational matrix C^AB has the components 

/ 1 \ 

1 / r) (3 cos i/j (3sm'ilj 

— (^/r)/3 cosV' 1 /3cos4' 

\ —f3smip —(3cosip 1 

Since the transformation flA4p is given by an orthogonal matrix the transformation be 
tween the corresponding sets of adjoint basis vectors is the same. 
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A2 The connection coefficients 

The non-trivial connection coefficients are given by the following expressions: 



r..C = a^{Z + i{r- iZ) K,) , T.,. = 



p p 



p p 



r = ^— (9 r/ — r 



p _ {rK2)' _ c K3 Q U 



r = ^— (9 w 



r — IL _ r 



(A6) 



r — iil _ £ ^2 



V?''' ii'j '^V' '''' 



where K^ = {Ks/K,){Ki/K3y. 

Other non-zero connection coefficients are obtained with help of antisymmetry over the first 
two indices: Tate = -^bac- 



APPENDIX B: PERTURBED EQUATIONS OF MOTION 

As we discuss in Section [3] our subset of equations describing the disc twist follows from 
equations (120|) after the procedure outlined in this Section. 
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The r-component of (!20|) gives 
K^m^p,J2U^ - ar^^ ^^^^^ + ^^^^ - ' 'rrr 



Ki r r 



dripU^v') + d^pU^v^ + ^""^JSJ pU^v' + F: = rd^p{U^fK^U + ^pUWd^W, (Bl) 

rA^il| ill 



where 

and the r, and ^-components give, respectively, 



Ki r 



K' Ki 
2^^ + a- ^ 



+ -f: 
p 



P 



K2 fiU^f + iU^f 



2ar-^U^ -^pv'P + ^^d^p, + 2—pd^v^+ 



dr{pU'^v') + d^pU^v^ + 



r^KiKl)' 
r^KiKl 



U^pv^ - a 



rKo 



-K2 . , e 



U^pv' + 



K-^rd^pU^U'^U + -piU'^fd^W, (B3) 
Kl r 



and 



K2 r.^^jr^r. . d^pi (U^y f KsU^\pi 

K——ru + U^o^v^ + r— h t 1 — 2ar———— h 

Ki ^ p r \ K2U'^ J p 

'K'l K'^ ar Kj 



2imv'^ 
'K2 



Kl 



Kl K2 2 K1K2 \K 



-d,U - 2a^ 



e KlZ [K 



+ -F^ 

p 

>.K, 
Kl 



rU^U^ + ar^^^U'^fd^U. (B4) 



.2 r K1K2 \K,, 

where are determined by the presence of viscous interactions in the disc. Since these 
interactions lead to a non-zero value of the drift component, U"^, we include terms containing 

in excepting the terms proportional to U^' through the dependence of on U^, see 
equation ( 14^ to keep the notation uniform. 

Explicitly, we have from (!20|) 

F: = —{d^Tf - rWd^T:^ - rd^pU^U^W, 



= {d^Tf - rWd^Tl'P) - rd^pW^WW, F; = S^TJ^, (B5) 



42 V. V. Zhuravlev and P. B. Ivanov 
and 

d.WiTl- + TJJ + a-§- d,Z {^T:- + ^T:x) , (B6) 

where = pW^U^. T^^ are components of the viscous stress tensor given by equations 
f l2T|) and f l22|) relevant for our purposes: 

T.^^ = (d^v^ - 2a^U^{U^fZ^ , = -r,r , (B7) 

where we use equations ( 1211) and ( l22l) . 

Note that we take into account all explicit terms of zero and first order in a in equations 
(IBltlB7l) . As we discuss in the main text, in fact, the only term linear in a, which is important 
for our purposes, is the gravitomagnetic term - the second term in the square brackets on 
the right hand side of flB4p . see equation flB12p below. Other terms may be important in 
other studies and are shown for a reference. The background quantities f/"^, U"^ and d^^p/ p 
can also be easily developed up to the linear in a order, and, for completeness, we show them 
below. We have 

where Ug and are given by (l25l) . and 

r n (f/v')2(]2 



p 



r2 ^2 



e, (B9) 



where 

n^ = Q- Qlt, (BIO) 

we remind that Vl = and the frequency of the Lense-Thirring precession (see e.g. Aliev 

& Galtsov 1981, Kato 1990), ^Ilt, is given by the expression 

^LT = = 2ai?~^ (Bll) 

It is important to note that from equation flB4p it follows that this characteristic fre- 
quency determines the disc's evolution in a certain regime. In order to show it let us tem- 
porarily assume that we consider a pressureless, inviscid disc with gas particles moving on 
circular orbits around a black hole. In this case we can set pi = pi = f ''^ = in flB4p . The 
last two terms on the right hand side can be also neglected since they contain products of 
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two small parameters: a{h/rY and a{td/ttw), respectively. Accordingly, we obtain from (IB4I1 

d^U = ^ltZ. (B12) 

Remembering that U = Z we see that in this case the disc rings precess with the precession 
frequency equal to Qlt- 



APPENDIX C: THE LAW OF CONSERVATION OF ANGULAR 
MOMENTUM 

The space-time of a Schwarzschild black hole is symmetric with respect to transformations 
belonging to the three-dimensional group of rotations, 0(3). Accordingly, equations of mo- 
tion contain three conserved quantities, which may be identified with three components of 
angular momentum. Conservation of the components perpendicular to the axis z requires 
a special divergent form of equation ( l48l) (or ( !6T|) ). In order to see that this equation has, 
indeed, such a form let us represent the term on the right hand side as —-pri^-§:;.F and the 
quantity F as F = Ci cos + sin (p, where Ci and C2 are functions of r and r. Introducing 
the complex notation used above we rewrite fHHj) as 



TV" 1 O 

Sf/^f/^(W - i^LT^) + -i(St/'^f/^ + f7)W' = F, (CI) 

K2 r^iil or 

where F = Ci + iC2- We multiply (IC1|) by r'^K'^ and integrate over r from some radius rj„ 
to rout > rin. We get 

L-?T + I = F(r,„)-F(r,„i), (C2) 

where 

ut p^out 

drr^K^m^U^W, T = / drr'^K^m'^U^^LT^ , (C3) 



and 

1= drr'^KiKl{m'PU'' + T'/)'W. (C4) 

It can be shown that L is proportional to a linear combination of x and y-components of 
the disc angular momentum within the region r^ < r < rout-, and T is the corresponding 
torque term provided by the gravitomagnetic force. 

The term I can be transformed with help of equation reflecting the law of conservation 
of angular momentum for a background disc model 

^ 'r^KiKl{m'<'U' + f'J)) + 2qU'P = 0, (C5) 



"'^KiKl dr 
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where q is flux of radiation from the disc surface o We integrate (1C4I) by parts and use ( ICSIl 
to get 

prout 

t-iT + 2 drr^K.KlqU^W = F(r,„) - F(r„„i) + I(r,„) - l(r,„0- (C6) 

Equation ( 1C6I) shows that when the black hole is not rotating (T = 0) a change of the 
angular momentum content within a region of the disc is determined by fluxes of angular 
momentum through the region inner and outer boundaries and angular momentum flux 
carried by radiation emitted from the disc surface. Thus, equation fHHj) has the required 
divergent form. 



This equation can either be obtained from equation l|20p written for a fiat disc or can be derived from the results of 
Page fc Thornj |l97J) and RH. 



